#### setting environment ####
require(foreign)
require(optiscale)
require(circular)
require(MNP)
require(coda)
require(ebal)
value.names <- c("Freedom", "Equality", "Economic stability", "Morality", 
                 "Self-reliance", "Social order", "Patriotism")
value.names.2 <- c("Freedom", "Equality", "Economic\nstability", "Morality", 
                   "Self-\nreliance", "Social\norder", "Patriotism")

# function to locate labels in Figure A1
label.location <- function(x) {
  loc <- rep(NA, length(x))
  for (i in 1:length(x)) {
    loc[i] <- sum(x[1:i]) - x[i] / 2
  }
  loc
}
# function to estimate the MDPREF model using the ALSOS algorithm
# This is nearly the same as a function in Jacoby (2014b).
ALSOS.MDPREF <- function(data, tolerance, max.iter) {
  vals.os <- data
  prev.fit <- 0
  niter <- 0
  improve <- 1
  while (improve > tolerance & niter <= max.iter) {
    niter <- niter + 1
    decomp <- svd(vals.os)
    d.sqd <- decomp$d ^ 2
    fit <- sum(d.sqd[1:2]) / sum(d.sqd)
    fitvector <- cumsum(d.sqd) / rep(sum(d.sqd), length(d.sqd))
    improve <- fit - prev.fit
    if (niter == 1) {
      history <- c(niter, fit, improve)
    } else {
      history <- rbind(history, c(niter, fit, improve))
    }
    subj.c1 <- decomp$u[, 1:2] %*% diag(decomp$d[1:2])
    pred.values <- subj.c1 %*% t(decomp$v[, 1:2])
    sumsqd <- apply(subj.c1 ^ 2, 1, sum)
    root.sumsqd <- (matrix(sumsqd, nrow = length(sumsqd), ncol = 1)) ^ 0.5
    subj.coord <- subj.c1 / (root.sumsqd %*% matrix(1, nrow = 1, ncol = 2))
    for (i in 1:nrow(vals.os)) {
      opped <- opscale(x.qual = data[i, ], x.quant = pred.values[i, ],
                       level = 2, process = 1, rescale = TRUE)
      vals.os[i,] <- opped$os
    }
    prev.fit <- fit
    print(paste0("iteration ", niter, " finished at ", date()))
  }
  result <- list(decomp = decomp, subj.coord = subj.coord, 
                 history = history, fitvector = fitvector)
  return(result)
}

#### preparing the datasets ####
## loading the UTAVS data
UTAVS <- read.csv("2014_2016UTASV20161004.csv")
UTAVS_DID <- read.csv("DID_population_ratios.csv")
# dealing with missing values
id <- UTAVS$ID
Q6 <- UTAVS$W1Q6
feeling <- UTAVS[, 23:33]
feeling2 <- UTAVS[, 112:119]
UTAVS[UTAVS == 66 | UTAVS == 99] <- NA
Q6[Q6 == 999] <- NA
feeling[feeling == 999] <- NA
feeling2[feeling2 == 999] <- NA
UTAVS$ID <- id
UTAVS$W1Q6 <- Q6
UTAVS[, 23:33] <- feeling
UTAVS[, 112:119] <- feeling2
# creating a variable for long-term partisanship
UTAVS <- data.frame(UTAVS, PID = ifelse(is.na(UTAVS$W1Q18_1), 11, UTAVS$W1Q18_1))
# add a variable for the municipal DID population ratios
UTAVS <- merge(UTAVS, UTAVS_DID, by = "ID", all.x = TRUE)

## loading the UTAES data
UTAES <- read.csv("2016UTASP20161010.csv")
# dealing with missing values
UTAES <- subset(UTAES, ! RESULT == 66)
district <- UTAES$DISTRICT
age <- UTAES$AGE
feeling <- UTAES[, 15:24]
UTAES[UTAES == 66 | UTAES == 99] <- NA
feeling[feeling == 999] <- NA
UTAES$DISTRICT <- district
UTAES$AGE <- age
UTAES[, 15:24] <- feeling
# creating a variable for party affiliation
party <- ifelse(UTAES$PARTY < 9, UTAES$PARTY, 
                ifelse(UTAES$PARTY == 11, 9, 
                       ifelse(UTAES$PARTY == 14, 10, 
                              ifelse(UTAES$PARTY == 15, 12, 11))))
UTAES <- data.frame(UTAES, party)
# data containing election winners only
UTAES.winners <- subset(UTAES, RESULT == 1)

## only use the data of those who completely ranked the seven values
UTAVS.complete <- subset(UTAVS, 
                         apply(UTAVS[, 67:73], 1, 
                               function(x) sum(x == 1) == 1 & sum(x == 2) == 1 & 
                                 sum(x == 3) == 1 & sum(x == 4) == 1 & 
                                 sum(x == 5) == 1 & sum(x == 6) == 1 & 
                                 sum(x == 7) == 1))
UTAVS.raw.values <- UTAVS.complete[, 67:73]
UTAVS.values <- 7 - UTAVS.raw.values
UTAVS.values.std <- t(apply(UTAVS.values, 1, scale))  # standadized value ranking for the MDPREF

UTAES.complete <- subset(UTAES, 
                         apply(UTAES[, 82:88], 1, 
                               function(x) sum(x == 1) == 1 & sum(x == 2) == 1 & 
                                 sum(x == 3) == 1 & sum(x == 4) == 1 & 
                                 sum(x == 5) == 1 & sum(x == 6) == 1 & 
                                 sum(x == 7) == 1))
UTAES.raw.values <- UTAES.complete[, 82:88]
UTAES.values <- 7 - UTAES.raw.values
UTAES.values.std <- t(apply(UTAES.values, 1, scale))  # standadized value ranking for the MDPREF

UTAES.complete.winners <- subset(UTAES.winners, 
                                 apply(UTAES.winners[, 82:88], 1, 
                                       function(x) sum(x == 1) == 1 & sum(x == 2) == 1 & 
                                         sum(x == 3) == 1 & sum(x == 4) == 1 & 
                                         sum(x == 5) == 1 & sum(x == 6) == 1 & 
                                         sum(x == 7) == 1))
UTAES.winners.raw.values <- UTAES.complete.winners[, 82:88]

## number of observations
nrow(UTAVS.complete)
nrow(UTAES.complete)
nrow(UTAES.complete.winners)

## Table A1 (patterns of missing responses by party in the UTAES)
missing.patterns <- matrix("", 13, 6)
missing.patterns[1:12, 1] <- sprintf("%3d", table(factor(UTAES$party, levels = c(1:12))))
missing.patterns[1:12, 2] <- sprintf("%3d", table(UTAES.complete$party))
missing.patterns[1:12, 3] <- sprintf("%3.0f", table(UTAES.complete$party) * 100 / table(UTAES$party))
missing.patterns[1:12, 4] <- sprintf("%3d", table(factor(UTAES.winners$party, levels = c(1:12))))
missing.patterns[1:12, 5] <- sprintf("%3d", table(factor(UTAES.complete.winners$party, levels = c(1:12))))
missing.patterns[1:12, 6] <- sprintf("%3.0f", table(factor(UTAES.complete.winners$party, levels = c(1:12))) * 100 / 
                                   table(factor(UTAES.winners$party, levels = c(1:12))), 1)
missing.patterns[13, 1] <- sprintf("%3d", nrow(UTAES))
missing.patterns[13, 2] <- sprintf("%3d", nrow(UTAES.complete))
missing.patterns[13, 3] <- sprintf("%3.0f", nrow(UTAES.complete) * 100 / nrow(UTAES))
missing.patterns[13, 4] <- sprintf("%3d", nrow(UTAES.winners))
missing.patterns[13, 5] <- sprintf("%3d", nrow(UTAES.complete.winners))
missing.patterns[13, 6] <- sprintf("%3.0f", nrow(UTAES.complete.winners) * 100 / nrow(UTAES.winners))
rownames(missing.patterns) <- c("LDP", "DP", "CGP", "JCP", "IfO", "SDP", "PLP&TYF", "JPK", "HRP", 
                                "Joint candidates of opposition parties", 
                                "Minor parties", "Other independents", "Total")
print(missing.patterns, quote = FALSE)

## number of inappropriate answers
# respondents who gave a partial ranking or did not rank any values
UTAVS.values.with.missing <- UTAVS[, 67:73]
sum(apply(UTAVS.values.with.missing, 1, function(x) sum(is.na(x))) > 0)
round(sum(apply(UTAVS.values.with.missing, 1, function(x) sum(is.na(x))) > 0) / nrow(UTAVS), 3)

UTAES.values.with.missing <- UTAES[, 82:88]
sum(apply(UTAES.values.with.missing, 1, function(x) sum(is.na(x))) > 0)
round(sum(apply(UTAES.values.with.missing, 1, function(x) sum(is.na(x))) > 0) / nrow(UTAES), 3)

#  respondents who ranked all seven values but gave a tied or non-sequential ranking
UTAVS.values.wo.missing <- na.omit(UTAVS[, 67:73])
nrow(UTAVS.values.wo.missing) - sum(apply(UTAVS.values.wo.missing, 1, 
                                          function(x) sum(x == 1) == 1 & sum(x == 2) == 1 & 
                                            sum(x == 3) == 1 & sum(x == 4) == 1 & 
                                            sum(x == 5) == 1 & sum(x == 6) == 1 & 
                                            sum(x == 7) == 1))
round((nrow(UTAVS.values.wo.missing) - sum(apply(UTAVS.values.wo.missing, 1, 
                                                 function(x) sum(x == 1) == 1 & sum(x == 2) == 1 & 
                                                   sum(x == 3) == 1 & sum(x == 4) == 1 & 
                                                   sum(x == 5) == 1 & sum(x == 6) == 1 & 
                                                   sum(x == 7) == 1))) / nrow(UTAVS), 3)

UTAES.values.wo.missing <- na.omit(UTAES[, 82:88])
nrow(UTAES.values.wo.missing) - sum(apply(UTAES.values.wo.missing, 1, 
                                          function(x) sum(x == 1) == 1 & sum(x == 2) == 1 & 
                                            sum(x == 3) == 1 & sum(x == 4) == 1 & 
                                            sum(x == 5) == 1 & sum(x == 6) == 1 & 
                                            sum(x == 7) == 1))
round((nrow(UTAES.values.wo.missing) - sum(apply(UTAES.values.wo.missing, 1, 
                                                 function(x) sum(x == 1) == 1 & sum(x == 2) == 1 & 
                                                   sum(x == 3) == 1 & sum(x == 4) == 1 & 
                                                   sum(x == 5) == 1 & sum(x == 6) == 1 & 
                                                   sum(x == 7) == 1))) / nrow(UTAES), 3)

# respondents who ranked the seven values in the same order as they were listed in the questionnaire 
sum(rowSums(abs(UTAVS.raw.values - tcrossprod(rep(1, nrow(UTAVS.raw.values)), c(1:7)))) == 0)
round(sum(rowSums(abs(UTAVS.raw.values - tcrossprod(rep(1, nrow(UTAVS.raw.values)), c(1:7)))) == 0) / nrow(UTAVS.raw.values), 3)

sum(rowSums(abs(UTAES.raw.values - tcrossprod(rep(1, nrow(UTAES.raw.values)), c(1:7)))) == 0)
round(sum(rowSums(abs(UTAES.raw.values - tcrossprod(rep(1, nrow(UTAES.raw.values)), c(1:7)))) == 0) / nrow(UTAES.raw.values), 3)

#### average ranking of values ####
## Figure 1 (average value ranks and the proportion choosing each value as their most important)
ranking.mean.UTAVS <- colMeans(UTAVS.raw.values)
no1.share.UTAVS <- apply(UTAVS.raw.values, 2, function(x) mean(x == 1)) * 100
round(apply(UTAVS.raw.values, 2, function(x) mean(x == 7)) * 100, 1)  # Which is the least prioritized value?

ranking.mean.UTAES <- colMeans(UTAES.raw.values)
no1.share.UTAES <- apply(UTAES.raw.values, 2, function(x) mean(x == 1)) * 100

ranking.mean.winners <- colMeans(UTAES.winners.raw.values)
no1.share.winners <- apply(UTAES.winners.raw.values, 2, function(x) mean(x == 1)) * 100

png("Figure_1.png", width = 5.5, height = 4.2, units = "in", pointsize = 8, res = 1500)
layout(matrix(c(1:4), 2, 2, byrow = TRUE))
par(mar=c(4, 4, 3, 1))
plot(ranking.mean.UTAVS, no1.share.UTAVS, xlim = c(2, 6), ylim = c(0, 55), 
     main = "Japanese Voters (2014)", 
     xlab = "Average Ranks", ylab = "Chosen as Most Important (%)")
text(ranking.mean.UTAVS[1], no1.share.UTAVS[1], "Freedom", pos = 3)
text(ranking.mean.UTAVS[2], no1.share.UTAVS[2], "Equality", pos = 3)
text(ranking.mean.UTAVS[3], no1.share.UTAVS[3], "Economic\nstability", pos = 3)
text(ranking.mean.UTAVS[4], no1.share.UTAVS[4], "Morality", pos = 3)
text(ranking.mean.UTAVS[5], no1.share.UTAVS[5], "Self-reliance", pos = 3)
text(ranking.mean.UTAVS[6], no1.share.UTAVS[6], "Social\norder", pos = 3)
text(ranking.mean.UTAVS[7] - 0.1, no1.share.UTAVS[7], "Patriotism", pos = 1)
plot(NULL, NULL, type = "n", xlim = c(2, 6), ylim = c(0, 55), 
     xlab = "", ylab = "", axes = FALSE)
plot(ranking.mean.UTAES, no1.share.UTAES, xlim = c(2, 6), ylim = c(0, 55), 
     main = "Japanese Candidates (2016)", 
     xlab = "Average Ranks", ylab = "Chosen as Most Important (%)")
text(ranking.mean.UTAES[1], no1.share.UTAES[1], "Freedom", pos = 1)
text(ranking.mean.UTAES[2], no1.share.UTAES[2], "Equality", pos = 1)
text(ranking.mean.UTAES[3], no1.share.UTAES[3], "Economic\nstability", pos = 3)
text(ranking.mean.UTAES[4], no1.share.UTAES[4], "Morality", pos = 4)
text(ranking.mean.UTAES[5], no1.share.UTAES[5], "Self-reliance", pos = 2)
text(ranking.mean.UTAES[6], no1.share.UTAES[6], "Social order", pos = 4)
text(ranking.mean.UTAES[7], no1.share.UTAES[7], "Patriotism", pos = 3)
plot(ranking.mean.winners, no1.share.winners, xlim = c(2, 6), ylim = c(0, 55), 
     main = "Japanese Members of the HoC (2016)", 
     xlab = "Average Ranks", ylab = "Chosen as Most Important (%)")
text(ranking.mean.winners[1], no1.share.winners[1], "Freedom", pos = 3)
text(ranking.mean.winners[2], no1.share.winners[2], "Equality", pos = 3)
text(ranking.mean.winners[3], no1.share.winners[3], "Economic\nstability", pos = 3)
text(ranking.mean.winners[4], no1.share.winners[4], "Morality", pos = 3)
text(ranking.mean.winners[5], no1.share.winners[5], "Self-reliance", pos = 2)
text(ranking.mean.winners[6], no1.share.winners[6], "Social\norder", pos = 1)
text(ranking.mean.winners[7] - 0.3, no1.share.winners[7], "Patriotism", pos = 3)
dev.off()

## Figure A1 (distribution of the ranks of values)
UTAVS.distribution <- apply(UTAVS.raw.values[, 7:1], 2, function(x) prop.table(table(x)))
UTAES.distribution <- apply(UTAES.raw.values[, 7:1], 2, function(x) prop.table(table(x)))
winners.distribution <- apply(UTAES.winners.raw.values[, 7:1], 2, function(x) prop.table(table(x)))
colnames(UTAVS.distribution) <- colnames(UTAES.distribution) <- 
  colnames(winners.distribution) <- c("Pat", "Soc", "Rel", "Mor", "Eco", "Equ", "Fre")

png("Figure_A1.png", width = 5.5, height = 5.5, units = "in", pointsize = 8, res = 1500)
layout(matrix(c(1:4), 2, 2, byrow = TRUE))
par(mar=c(4, 3, 3, 1))
barplot(UTAVS.distribution, horiz = TRUE, main = "Japanese Voters (2014)", 
        axes = FALSE, cex.names = 0.9)
axis(1, at = c(0, 0.25, 0.5, 0.75, 1), labels = c("0", "25", "50", "75", "100"))
mtext("Cumulative Proportion (%)", side = 1, line = 2.5)
mtext(c(1:7), line = -0.5, at = label.location(UTAVS.distribution[, 7]), cex = 0.8)
plot(NULL, NULL, type = "n", xlim = c(0, 1), ylim = c(0, 1), 
     xlab = "", ylab = "", axes = FALSE)
barplot(UTAES.distribution, horiz = TRUE, main = "Japanese Candidates (2016)", 
        axes = FALSE, cex.names = 0.9)
axis(1, at = c(0, 0.25, 0.5, 0.75, 1), labels = c("0", "25", "50", "75", "100"))
mtext("Cumulative Proportion (%)", side = 1, line = 2.5)
mtext(c(1:7), line = -0.5, at = label.location(UTAES.distribution[, 7]), cex = 0.8)
barplot(winners.distribution, horiz = TRUE, main = "Japanese Members of the HoC (2016)", 
        axes = FALSE, cex.names = 0.9)
axis(1, at = c(0, 0.25, 0.5, 0.75, 1), labels = c("0", "25", "50", "75", "100"))
mtext("Cumulative Proportion (%)", side = 1, line = 2.5)
mtext(c(1:7), line = -0.5, at = label.location(winners.distribution[, 7]), cex = 0.8)
dev.off()

#### value structures ####
## estimate the MDPREF model
# MDPREF fot the UTAVS
mdpref.UTAVS <- ALSOS.MDPREF(data = UTAVS.values.std, tolerance = 0.01, max.iter = 25)
# R squared values
round(mdpref.UTAVS$fitvector, 3)
# rotate the results
raw.freedom <- as.numeric(coord2rad(mdpref.UTAVS$decomp$v[1, 1], mdpref.UTAVS$decomp$v[1, 2]))
UTAVS.rotation.matrix <- matrix(c(cos(raw.freedom - pi / 2), sin(raw.freedom - pi / 2), 
                                  -sin(raw.freedom - pi / 2), cos(raw.freedom - pi / 2)), 2, 2)
UTAVS.value.coords <- mdpref.UTAVS$decomp$v[, 1:2] %*% UTAVS.rotation.matrix
rownames(UTAVS.value.coords) <- value.names
UTAVS.obs.coords <- mdpref.UTAVS$subj.coord %*% UTAVS.rotation.matrix

# MDPREF fot the UTAES
mdpref.UTAES <- ALSOS.MDPREF(data = UTAES.values.std, tolerance = 0.01, max.iter = 25)
# R squared values
round(mdpref.UTAES$fitvector, 3)
# rotate the results
raw.freedom <- as.numeric(coord2rad(mdpref.UTAES$decomp$v[1, 1], mdpref.UTAES$decomp$v[1, 2]))
UTAES.rotation.matrix <- matrix(c(cos(raw.freedom - pi / 2), sin(raw.freedom - pi / 2), 
                                  -sin(raw.freedom - pi / 2), cos(raw.freedom - pi / 2)), 2, 2)
UTAES.value.coords <- mdpref.UTAES$decomp$v[, 1:2] %*% UTAES.rotation.matrix
rownames(UTAES.value.coords) <- value.names
UTAES.obs.coords <- mdpref.UTAES$subj.coord %*% UTAES.rotation.matrix

## Figure 2 (illustration of an estimation result interpretation using MDPREF)
point.A <- c(cos(0.8 * pi), sin(0.8 * pi))
projection.A <- tcrossprod(UTAVS.value.coords[c(1, 6, 7), ] %*% point.A, point.A)
point.B <- c(cos(1.4 * pi), sin(1.4 * pi))
projection.B <- tcrossprod(UTAVS.value.coords[c(1, 6, 7), ] %*% point.B, point.B)

png("Figure_2.png", width = 4, height = 4, units = "in", pointsize = 9, res = 1500)
par(mar=c(0, 0, 0, 0))
plot(UTAVS.value.coords[c(1, 6, 7), 1], UTAVS.value.coords[c(1, 6, 7), 2], type = "p", bty = "n", 
     xlim = c(-1, 1), ylim = c(-1, 1), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
text(UTAVS.value.coords[c(1, 6, 7), 1], UTAVS.value.coords[c(1, 6, 7), 2], 
     labels = value.names.2[c(1, 6, 7)], pos = c(3, 2, 4))
points(0, 0, pch = 19)
text(0, 0, labels = "O", pos = 4)
points(point.A[1], point.A[2], pch = 19)
text(point.A[1], point.A[2], labels = "A", pos = 2)
arrows(-point.A[1] * 0.97, -point.A[2] * 0.97, point.A[1] * 0.97, point.A[2] * 0.97, length = 0.05)
segments(UTAVS.value.coords[c(1, 6, 7), 1], UTAVS.value.coords[c(1, 6, 7), 2], 
         projection.A[, 1], projection.A[, 2], lty = 3)
text(projection.A[, 1], projection.A[, 2], labels = c(1:3), cex = 0.8, pos = c(1, 3, 1))
points(point.B[1], point.B[2], pch = 19)
text(point.B[1], point.B[2], labels = "B", pos = 1)
arrows(-point.B[1] * 0.97, -point.B[2] * 0.97, point.B[1] * 0.97, point.B[2] * 0.97, length = 0.05)
segments(UTAVS.value.coords[c(1, 6, 7), 1], UTAVS.value.coords[c(1, 6, 7), 2], 
         projection.B[, 1], projection.B[, 2], lty = 3)
text(projection.B[, 1], projection.B[, 2], labels = c(3, 1, 2), cex = 0.8, pos = c(4, 4, 2))
dev.off()

## Figure 3 (value structures and the distribution of value preferences)
png("Figure_3.png", width = 6.4, height = 6.4, units = "in", pointsize = 11, res = 1500)
layout(matrix(c(1:4), 2, 2, byrow = TRUE))
par(mar=c(0.5, 0.5, 3, 0.5))
set.seed(12345)
plot(density.circular(coord2rad(UTAVS.obs.coords[, 1], UTAVS.obs.coords[, 2]), 
                      kernel = "vonmises", bw = 25), main = "Japanese Voters (2014)", axes = FALSE, shrink = 1.2)
points(UTAVS.value.coords[, 1], UTAVS.value.coords[, 2])
text(UTAVS.value.coords[, 1], UTAVS.value.coords[, 2], labels = value.names, cex = 0.8, pos = 3)
arrows.circular(coord2rad(mean(UTAVS.obs.coords[, 1]), mean(UTAVS.obs.coords[, 2])), 
                y = (mean(UTAVS.obs.coords[, 1]) ^ 2 + mean(UTAVS.obs.coords[, 2]) ^ 2) ^ 0.5, 
                length = 0.05)
sampled.id <- sample(c(1:nrow(UTAVS.complete)), 300)
ticks.circular(coord2rad(UTAVS.obs.coords[, 1], UTAVS.obs.coords[, 2])[sampled.id])
text(1.2, -1.2, labels = paste0("N = ", nrow(UTAVS.obs.coords)))
plot(NULL, NULL, type = "n", xlim = c(0, 1), ylim = c(0, 1), 
     xlab = "", ylab = "", axes = FALSE)
plot(density.circular(coord2rad(UTAES.obs.coords[, 1], UTAES.obs.coords[, 2]), 
                      kernel = "vonmises", bw = 25), main = "Japanese Candidates (2016)", axes = FALSE, shrink = 1.2)
points(UTAES.value.coords[, 1], UTAES.value.coords[, 2])
text(UTAES.value.coords[, 1], UTAES.value.coords[, 2], labels = value.names.2, cex = 0.8, pos = c(3, 2, 3, 3, 3, 3, 3))
arrows.circular(coord2rad(mean(UTAES.obs.coords[, 1]), mean(UTAES.obs.coords[, 2])), 
                y = (mean(UTAES.obs.coords[, 1]) ^ 2 + mean(UTAES.obs.coords[, 2]) ^ 2) ^ 0.5, 
                length = 0.05)
ticks.circular(coord2rad(UTAES.obs.coords[, 1], UTAES.obs.coords[, 2]))
text(1.2, -1.2, labels = paste0("N = ", nrow(UTAES.obs.coords)))
plot(density.circular(coord2rad(UTAES.obs.coords[UTAES.complete$RESULT == 1, 1], 
                                UTAES.obs.coords[UTAES.complete$RESULT == 1, 2]), 
                      kernel = "vonmises", bw = 25), main = "Japanese Members of the HoC (2016)", axes = FALSE, shrink = 1.2)
points(UTAES.value.coords[, 1], UTAES.value.coords[, 2])
text(UTAES.value.coords[, 1], UTAES.value.coords[, 2], labels = value.names.2, cex = 0.8, pos = c(3, 4, 3, 3, 3, 3, 3))
arrows.circular(coord2rad(mean(UTAES.obs.coords[UTAES.complete$RESULT == 1, 1]), 
                          mean(UTAES.obs.coords[UTAES.complete$RESULT == 1, 2])), 
                y = (mean(UTAES.obs.coords[UTAES.complete$RESULT == 1, 1]) ^ 2 + 
                       mean(UTAES.obs.coords[UTAES.complete$RESULT == 1, 2]) ^ 2) ^ 0.5, 
                length = 0.05)
ticks.circular(coord2rad(UTAES.obs.coords[UTAES.complete$RESULT == 1, 1], 
                         UTAES.obs.coords[UTAES.complete$RESULT == 1, 2]))
text(1.2, -1.2, labels = paste0("N = ", sum(UTAES.complete$RESULT == 1)))
dev.off()

#### party-based comparison among Japanese voters (overall value preferences) ####
## Figure 4 (value structures and the distribution of value preferences of Japanese voters according to long-term partisanship)
party.labels.UTAVS <- c("LDP", "DPJ", "JIP", "CGP", "PFG", "JCP", "PLP", "SDP", "", "", "Independent")
png("Figure_4.png", width = 4.8, height = 4.8, units = "in", pointsize = 9, res = 1500)
layout(matrix(c(1:9), 3, 3, byrow = TRUE))
par(mar=c(0, 0, 2, 0))
for (i in c(1:8, 11)) {
  plot(density.circular(coord2rad(UTAVS.obs.coords[UTAVS.complete$PID == i, 1], 
                                  UTAVS.obs.coords[UTAVS.complete$PID == i, 2]), 
                        kernel = "vonmises", bw = 25), 
       main = party.labels.UTAVS[i], axes = FALSE, shrink = 1.3, 
       xlim = c(-1.3, 1.3), ylim = c(-1.15, 1.45))
  points(UTAVS.value.coords[, 1], UTAVS.value.coords[, 2])
  points(0, 0, pch = 20, cex = 0.5)
  arrows.circular(coord2rad(mean(UTAVS.obs.coords[UTAVS.complete$PID == i, 1]), 
                            mean(UTAVS.obs.coords[UTAVS.complete$PID == i, 2])), 
                  y = (mean(UTAVS.obs.coords[UTAVS.complete$PID == i, 1]) ^ 2 + 
                         mean(UTAVS.obs.coords[UTAVS.complete$PID == i, 2]) ^ 2) ^ 0.5, 
                  length = 0.03)
  text(1.3, -1.15, labels = paste0("N = ", sum(UTAVS.complete$PID == i)))
}
dev.off()

#### party-based comparison among Japanese voters (rank order of individual values) ####
## preparing the data and variables
# dependent variable (inconsistent ties and subsequent values are corrected)
all.values <- matrix(NA, nrow(UTAVS), 7)
for (i in 1:nrow(UTAVS)) {
  count <- 1
  ranking <- UTAVS[i, c(67:73)]
  if (sum(ranking == 7, na.rm = TRUE) > 0) {
    all.values[i, which(ranking == 7)] <- count
    count <- count + 1
  }
  if (sum(ranking == 6, na.rm = TRUE) > 0) {
    all.values[i, which(ranking == 6)] <- count
    count <- count + 1
  }
  if (sum(ranking == 5, na.rm = TRUE) > 0) {
    all.values[i, which(ranking == 5)] <- count
    count <- count + 1
  }
  if (sum(ranking == 4, na.rm = TRUE) > 0) {
    all.values[i, which(ranking == 4)] <- count
    count <- count + 1
  }
  if (sum(ranking == 3, na.rm = TRUE) > 0) {
    all.values[i, which(ranking == 3)] <- count
    count <- count + 1
  }
  if (sum(ranking == 2, na.rm = TRUE) > 0) {
    all.values[i, which(ranking == 2)] <- count
    count <- count + 1
  }
  if (sum(ranking == 1, na.rm = TRUE) > 0) {
    all.values[i, which(ranking == 1)] <- count
    count <- count + 1
  }
}
colnames(all.values) <- c("freedom", "equality", "econ.sta", "morality", 
                          "self.reli", "soc.order", "patriotism")

# independent variables
LDP <- ifelse(UTAVS$PID == 1, 1, 0)
DPJ <- ifelse(UTAVS$PID == 2, 1, 0)
JRP <- ifelse(UTAVS$PID == 3, 1, 0)
female <- UTAVS$W1F1 - 1
age <- UTAVS$W1F2
education <- ifelse(UTAVS$W1F3 == 1, 1, 
                    ifelse(UTAVS$W1F3 == 2, 2, 
                           ifelse(UTAVS$W1F3 < 5, 3, 
                                  ifelse(UTAVS$W1F3 < 7, 4, NA))))
DID <- UTAVS$DID_CITY

# dataset
number.of.NA <- apply(all.values, 1, function(x) sum(is.na(x)))
MVNOS.data <- data.frame(all.values, female = female, age = age, 
                       education = education, DID = DID, 
                       LDP, DPJ, JRP)[number.of.NA < 6, ]  # respondents who provided rank orders for at least two values
nrow(MVNOS.data)  # number of observations

## estimate the MVNOS model without demographic variables
MVNOS.wo.demographics <- list()
set.seed(123)
MVNOS.wo.demographics[[1]] <- mnp(cbind(freedom, equality, econ.sta, morality, 
                                      self.reli, soc.order, patriotism) ~ LDP + DPJ + JRP, 
                                data = MVNOS.data, n.draws = 11000, burnin = 1000, thin = 9)
set.seed(456)
MVNOS.wo.demographics[[2]] <- mnp(cbind(freedom, equality, econ.sta, morality, 
                                      self.reli, soc.order, patriotism) ~ LDP + DPJ + JRP, 
                                data = MVNOS.data, n.draws = 11000, burnin = 1000, thin = 9)
set.seed(789)
MVNOS.wo.demographics[[3]] <- mnp(cbind(freedom, equality, econ.sta, morality, 
                                      self.reli, soc.order, patriotism) ~ LDP + DPJ + JRP, 
                                data = MVNOS.data, n.draws = 11000, burnin = 1000, thin = 9)
#save(MVNOS.wo.demographics, file = "MVNOS_without_demographics.Rdata")
#load("MVNOS_without_demographics.Rdata")

# convergence check
which(gelman.diag(mcmc.list(as.mcmc(MVNOS.wo.demographics[[1]]$param), 
                            as.mcmc(MVNOS.wo.demographics[[2]]$param), 
                            as.mcmc(MVNOS.wo.demographics[[3]]$param)), 
                  multivariate = FALSE)[[1]][, 1] > 1.1)

## post-estimation simulations without demographic variables
# datasets used for simulations
newdata.LDP <- newdata.DPJ <- newdata.JRP <- MVNOS.data
newdata.LDP[, 12:14] <- newdata.DPJ[, 12:14] <- newdata.JRP[, 12:14] <- 0
newdata.LDP$LDP <- 1
newdata.DPJ$DPJ <- 1
newdata.JRP$JRP <- 1

# prediction, manipulating voters' long-term partisanship
pred.MVNOS.LDP <- pred.MVNOS.DPJ <- pred.MVNOS.JRP <- list()
set.seed(12345)
pred.MVNOS.LDP[[1]] <- predict(MVNOS.wo.demographics[[1]], newdata = newdata.LDP, type = "order")
pred.MVNOS.LDP[[2]] <- predict(MVNOS.wo.demographics[[2]], newdata = newdata.LDP, type = "order")
pred.MVNOS.LDP[[3]] <- predict(MVNOS.wo.demographics[[3]], newdata = newdata.LDP, type = "order")
pred.MVNOS.DPJ[[1]] <- predict(MVNOS.wo.demographics[[1]], newdata = newdata.DPJ, type = "order")
pred.MVNOS.DPJ[[2]] <- predict(MVNOS.wo.demographics[[2]], newdata = newdata.DPJ, type = "order")
pred.MVNOS.DPJ[[3]] <- predict(MVNOS.wo.demographics[[3]], newdata = newdata.DPJ, type = "order")
pred.MVNOS.JRP[[1]] <- predict(MVNOS.wo.demographics[[1]], newdata = newdata.JRP, type = "order")
pred.MVNOS.JRP[[2]] <- predict(MVNOS.wo.demographics[[2]], newdata = newdata.JRP, type = "order")
pred.MVNOS.JRP[[3]] <- predict(MVNOS.wo.demographics[[3]], newdata = newdata.JRP, type = "order")

# compute differences in the average rank orders
dif.MVNOS.DPJ.LDP <- cbind(colMeans(rbind(t(apply(pred.MVNOS.DPJ[[1]]$o, c(2, 3), mean)), 
                                        t(apply(pred.MVNOS.DPJ[[2]]$o, c(2, 3), mean)), 
                                        t(apply(pred.MVNOS.DPJ[[3]]$o, c(2, 3), mean))) - 
                                    rbind(t(apply(pred.MVNOS.LDP[[1]]$o, c(2, 3), mean)), 
                                          t(apply(pred.MVNOS.LDP[[2]]$o, c(2, 3), mean)), 
                                          t(apply(pred.MVNOS.LDP[[3]]$o, c(2, 3), mean)))), 
                         HPDinterval(as.mcmc(rbind(t(apply(pred.MVNOS.DPJ[[1]]$o, c(2, 3), mean)), 
                                                   t(apply(pred.MVNOS.DPJ[[2]]$o, c(2, 3), mean)), 
                                                   t(apply(pred.MVNOS.DPJ[[3]]$o, c(2, 3), mean))) - 
                                               rbind(t(apply(pred.MVNOS.LDP[[1]]$o, c(2, 3), mean)), 
                                                     t(apply(pred.MVNOS.LDP[[2]]$o, c(2, 3), mean)), 
                                                     t(apply(pred.MVNOS.LDP[[3]]$o, c(2, 3), mean))))))
dif.MVNOS.JRP.LDP <- cbind(colMeans(rbind(t(apply(pred.MVNOS.JRP[[1]]$o, c(2, 3), mean)), 
                                        t(apply(pred.MVNOS.JRP[[2]]$o, c(2, 3), mean)), 
                                        t(apply(pred.MVNOS.JRP[[3]]$o, c(2, 3), mean))) - 
                                    rbind(t(apply(pred.MVNOS.LDP[[1]]$o, c(2, 3), mean)), 
                                          t(apply(pred.MVNOS.LDP[[2]]$o, c(2, 3), mean)), 
                                          t(apply(pred.MVNOS.LDP[[3]]$o, c(2, 3), mean)))), 
                         HPDinterval(as.mcmc(rbind(t(apply(pred.MVNOS.JRP[[1]]$o, c(2, 3), mean)), 
                                                   t(apply(pred.MVNOS.JRP[[2]]$o, c(2, 3), mean)), 
                                                   t(apply(pred.MVNOS.JRP[[3]]$o, c(2, 3), mean))) - 
                                               rbind(t(apply(pred.MVNOS.LDP[[1]]$o, c(2, 3), mean)), 
                                                     t(apply(pred.MVNOS.LDP[[2]]$o, c(2, 3), mean)), 
                                                     t(apply(pred.MVNOS.LDP[[3]]$o, c(2, 3), mean))))))
dif.MVNOS.JRP.DPJ <- cbind(colMeans(rbind(t(apply(pred.MVNOS.JRP[[1]]$o, c(2, 3), mean)), 
                                        t(apply(pred.MVNOS.JRP[[2]]$o, c(2, 3), mean)), 
                                        t(apply(pred.MVNOS.JRP[[3]]$o, c(2, 3), mean))) - 
                                    rbind(t(apply(pred.MVNOS.DPJ[[1]]$o, c(2, 3), mean)), 
                                          t(apply(pred.MVNOS.DPJ[[2]]$o, c(2, 3), mean)), 
                                          t(apply(pred.MVNOS.DPJ[[3]]$o, c(2, 3), mean)))), 
                         HPDinterval(as.mcmc(rbind(t(apply(pred.MVNOS.JRP[[1]]$o, c(2, 3), mean)), 
                                                   t(apply(pred.MVNOS.JRP[[2]]$o, c(2, 3), mean)), 
                                                   t(apply(pred.MVNOS.JRP[[3]]$o, c(2, 3), mean))) - 
                                               rbind(t(apply(pred.MVNOS.DPJ[[1]]$o, c(2, 3), mean)), 
                                                     t(apply(pred.MVNOS.DPJ[[2]]$o, c(2, 3), mean)), 
                                                     t(apply(pred.MVNOS.DPJ[[3]]$o, c(2, 3), mean))))))

# 90% credible intervals
HPDinterval(as.mcmc(rbind(t(apply(pred.MVNOS.JRP[[1]]$o, c(2, 3), mean)), 
                          t(apply(pred.MVNOS.JRP[[2]]$o, c(2, 3), mean)), 
                          t(apply(pred.MVNOS.JRP[[3]]$o, c(2, 3), mean))) - 
                      rbind(t(apply(pred.MVNOS.LDP[[1]]$o, c(2, 3), mean)), 
                            t(apply(pred.MVNOS.LDP[[2]]$o, c(2, 3), mean)), 
                            t(apply(pred.MVNOS.LDP[[3]]$o, c(2, 3), mean)))), prob = 0.9)
HPDinterval(as.mcmc(rbind(t(apply(pred.MVNOS.JRP[[1]]$o, c(2, 3), mean)), 
                          t(apply(pred.MVNOS.JRP[[2]]$o, c(2, 3), mean)), 
                          t(apply(pred.MVNOS.JRP[[3]]$o, c(2, 3), mean))) - 
                      rbind(t(apply(pred.MVNOS.DPJ[[1]]$o, c(2, 3), mean)), 
                            t(apply(pred.MVNOS.DPJ[[2]]$o, c(2, 3), mean)), 
                            t(apply(pred.MVNOS.DPJ[[3]]$o, c(2, 3), mean)))), prob = 0.9)

## Figure 5 (simulation results on the relation between ranking of values and long-term partisanship)
png("Figure_5.png", width = 5, height = 3.4, units = "in", pointsize = 8, res = 1500)
layout(matrix(c(1:4), 2, 2, byrow = TRUE))
par(mar=c(4, 1, 3, 1))
dotchart(dif.MVNOS.DPJ.LDP[c(1, 7:2), 1], labels = rev(value.names), pch = 20, xlim = c(-1, 1), 
         main = "LDP \u2192 DPJ", xlab = "Change in Average Rank")
points(dif.MVNOS.DPJ.LDP[c(1, 7:2), 1], c(1:7), pch = 19)
segments(dif.MVNOS.DPJ.LDP[c(1, 7:2), 2], c(1:7), dif.MVNOS.DPJ.LDP[c(1, 7:2), 3], c(1:7))
abline(v = 0, lty = 3)
dotchart(dif.MVNOS.JRP.LDP[c(1, 7:2), 1], labels = rev(value.names), pch = 20, xlim = c(-1, 1), 
         main = "LDP \u2192 JIP", xlab = "Change in Average Rank")
points(dif.MVNOS.JRP.LDP[c(1, 7:2), 1], c(1:7), pch = 19)
segments(dif.MVNOS.JRP.LDP[c(1, 7:2), 2], c(1:7), dif.MVNOS.JRP.LDP[c(1, 7:2), 3], c(1:7))
abline(v = 0, lty = 3)
dotchart(dif.MVNOS.JRP.DPJ[c(1, 7:2), 1], labels = rev(value.names), pch = 20, xlim = c(-1, 1), 
         main = "DPJ \u2192 JIP", xlab = "Change in Average Rank")
points(dif.MVNOS.JRP.DPJ[c(1, 7:2), 1], c(1:7), pch = 19)
segments(dif.MVNOS.JRP.DPJ[c(1, 7:2), 2], c(1:7), dif.MVNOS.JRP.DPJ[c(1, 7:2), 3], c(1:7))
abline(v = 0, lty = 3)
dev.off()

## estimate the MVNOS model with demographic variables
MVNOS.with.demographics <- list()
set.seed(123)
MVNOS.with.demographics[[1]] <- mnp(cbind(freedom, equality, econ.sta, morality, 
                                        self.reli, soc.order, patriotism) ~ 
                                    female + age + education + DID + LDP + DPJ + JRP, 
                                  data = MVNOS.data, n.draws = 11000, burnin = 1000, thin = 9)
set.seed(456)
MVNOS.with.demographics[[2]] <- mnp(cbind(freedom, equality, econ.sta, morality, 
                                        self.reli, soc.order, patriotism) ~ 
                                    female + age + education + DID + LDP + DPJ + JRP, 
                                  data = MVNOS.data, n.draws = 11000, burnin = 1000, thin = 9)
set.seed(789)
MVNOS.with.demographics[[3]] <- mnp(cbind(freedom, equality, econ.sta, morality, 
                                        self.reli, soc.order, patriotism) ~ 
                                    female + age + education + DID + LDP + DPJ + JRP, 
                                  data = MVNOS.data, n.draws = 11000, burnin = 1000, thin = 9)
#save(MVNOS.with.demographics, file = "MVNOS_with_demographics.Rdata")
#load("MVNOS_with_demographics.Rdata")

# convergence check
which(gelman.diag(mcmc.list(as.mcmc(MVNOS.with.demographics[[1]]$param), 
                            as.mcmc(MVNOS.with.demographics[[2]]$param), 
                            as.mcmc(MVNOS.with.demographics[[3]]$param)), multivariate = FALSE)[[1]][, 1] > 1.1)

# number of observations after respondents with missing independent variables were discarded 
nrow(MVNOS.data)
nrow(MVNOS.with.demographics[[1]]$x)
nrow(MVNOS.data) - nrow(MVNOS.with.demographics[[1]]$x)
round((nrow(MVNOS.data) - nrow(MVNOS.with.demographics[[1]]$x)) / nrow(MVNOS.data), 3)

## post-estimation simulations with demographic variables
# female
newdata.male <- newdata.female <- MVNOS.data
newdata.male$female <- 0
newdata.female$female <- 1
pred.MVNOS.male <- pred.MVNOS.female <- list()
set.seed(1000)
pred.MVNOS.male[[1]] <- predict(MVNOS.with.demographics[[1]], newdata = newdata.male, type = "order")
pred.MVNOS.male[[2]] <- predict(MVNOS.with.demographics[[2]], newdata = newdata.male, type = "order")
pred.MVNOS.male[[3]] <- predict(MVNOS.with.demographics[[3]], newdata = newdata.male, type = "order")
pred.MVNOS.female[[1]] <- predict(MVNOS.with.demographics[[1]], newdata = newdata.female, type = "order")
pred.MVNOS.female[[2]] <- predict(MVNOS.with.demographics[[2]], newdata = newdata.female, type = "order")
pred.MVNOS.female[[3]] <- predict(MVNOS.with.demographics[[3]], newdata = newdata.female, type = "order")
dif.MVNOS.female <- cbind(colMeans(rbind(t(apply(pred.MVNOS.female[[1]]$o, c(2, 3), mean)), 
                                       t(apply(pred.MVNOS.female[[2]]$o, c(2, 3), mean)), 
                                       t(apply(pred.MVNOS.female[[3]]$o, c(2, 3), mean))) - 
                                   rbind(t(apply(pred.MVNOS.male[[1]]$o, c(2, 3), mean)), 
                                         t(apply(pred.MVNOS.male[[2]]$o, c(2, 3), mean)), 
                                         t(apply(pred.MVNOS.male[[3]]$o, c(2, 3), mean)))), 
                        HPDinterval(as.mcmc(rbind(t(apply(pred.MVNOS.female[[1]]$o, c(2, 3), mean)), 
                                                  t(apply(pred.MVNOS.female[[2]]$o, c(2, 3), mean)), 
                                                  t(apply(pred.MVNOS.female[[3]]$o, c(2, 3), mean))) - 
                                              rbind(t(apply(pred.MVNOS.male[[1]]$o, c(2, 3), mean)), 
                                                    t(apply(pred.MVNOS.male[[2]]$o, c(2, 3), mean)), 
                                                    t(apply(pred.MVNOS.male[[3]]$o, c(2, 3), mean))))))

# age
newdata.young <- newdata.old <- MVNOS.data
newdata.young$age <- 1
newdata.old$age <- 6
pred.MVNOS.young <- pred.MVNOS.old <- list()
set.seed(2000)
pred.MVNOS.young[[1]] <- predict(MVNOS.with.demographics[[1]], newdata = newdata.young, type = "order")
pred.MVNOS.young[[2]] <- predict(MVNOS.with.demographics[[2]], newdata = newdata.young, type = "order")
pred.MVNOS.young[[3]] <- predict(MVNOS.with.demographics[[3]], newdata = newdata.young, type = "order")
pred.MVNOS.old[[1]] <- predict(MVNOS.with.demographics[[1]], newdata = newdata.old, type = "order")
pred.MVNOS.old[[2]] <- predict(MVNOS.with.demographics[[2]], newdata = newdata.old, type = "order")
pred.MVNOS.old[[3]] <- predict(MVNOS.with.demographics[[3]], newdata = newdata.old, type = "order")
dif.MVNOS.age <- cbind(colMeans(rbind(t(apply(pred.MVNOS.old[[1]]$o, c(2, 3), mean)), 
                                    t(apply(pred.MVNOS.old[[2]]$o, c(2, 3), mean)), 
                                    t(apply(pred.MVNOS.old[[3]]$o, c(2, 3), mean))) - 
                                rbind(t(apply(pred.MVNOS.young[[1]]$o, c(2, 3), mean)), 
                                      t(apply(pred.MVNOS.young[[2]]$o, c(2, 3), mean)), 
                                      t(apply(pred.MVNOS.young[[3]]$o, c(2, 3), mean)))), 
                     HPDinterval(as.mcmc(rbind(t(apply(pred.MVNOS.old[[1]]$o, c(2, 3), mean)), 
                                               t(apply(pred.MVNOS.old[[2]]$o, c(2, 3), mean)), 
                                               t(apply(pred.MVNOS.old[[3]]$o, c(2, 3), mean))) - 
                                           rbind(t(apply(pred.MVNOS.young[[1]]$o, c(2, 3), mean)), 
                                                 t(apply(pred.MVNOS.young[[2]]$o, c(2, 3), mean)), 
                                                 t(apply(pred.MVNOS.young[[3]]$o, c(2, 3), mean))))))

# education
newdata.primary <- newdata.college <- MVNOS.data
newdata.primary$education <- 1
newdata.college$education <- 4
pred.MVNOS.primary <- pred.MVNOS.college <- list()
set.seed(3000)
pred.MVNOS.primary[[1]] <- predict(MVNOS.with.demographics[[1]], newdata = newdata.primary, type = "order")
pred.MVNOS.primary[[2]] <- predict(MVNOS.with.demographics[[2]], newdata = newdata.primary, type = "order")
pred.MVNOS.primary[[3]] <- predict(MVNOS.with.demographics[[3]], newdata = newdata.primary, type = "order")
pred.MVNOS.college[[1]] <- predict(MVNOS.with.demographics[[1]], newdata = newdata.college, type = "order")
pred.MVNOS.college[[2]] <- predict(MVNOS.with.demographics[[2]], newdata = newdata.college, type = "order")
pred.MVNOS.college[[3]] <- predict(MVNOS.with.demographics[[3]], newdata = newdata.college, type = "order")
dif.MVNOS.education <- cbind(colMeans(rbind(t(apply(pred.MVNOS.college[[1]]$o, c(2, 3), mean)), 
                                          t(apply(pred.MVNOS.college[[2]]$o, c(2, 3), mean)), 
                                          t(apply(pred.MVNOS.college[[3]]$o, c(2, 3), mean))) - 
                                      rbind(t(apply(pred.MVNOS.primary[[1]]$o, c(2, 3), mean)), 
                                            t(apply(pred.MVNOS.primary[[2]]$o, c(2, 3), mean)), 
                                            t(apply(pred.MVNOS.primary[[3]]$o, c(2, 3), mean)))), 
                           HPDinterval(as.mcmc(rbind(t(apply(pred.MVNOS.college[[1]]$o, c(2, 3), mean)), 
                                                     t(apply(pred.MVNOS.college[[2]]$o, c(2, 3), mean)), 
                                                     t(apply(pred.MVNOS.college[[3]]$o, c(2, 3), mean))) - 
                                                 rbind(t(apply(pred.MVNOS.primary[[1]]$o, c(2, 3), mean)), 
                                                       t(apply(pred.MVNOS.primary[[2]]$o, c(2, 3), mean)), 
                                                       t(apply(pred.MVNOS.primary[[3]]$o, c(2, 3), mean))))))

# DID population ratio
newdata.rural <- newdata.urban <- MVNOS.data
newdata.rural$DID <- 0
newdata.urban$DID <- 1
pred.MVNOS.rural <- pred.MVNOS.urban <- list()
set.seed(4000)
pred.MVNOS.rural[[1]] <- predict(MVNOS.with.demographics[[1]], newdata = newdata.rural, type = "order")
pred.MVNOS.rural[[2]] <- predict(MVNOS.with.demographics[[2]], newdata = newdata.rural, type = "order")
pred.MVNOS.rural[[3]] <- predict(MVNOS.with.demographics[[3]], newdata = newdata.rural, type = "order")
pred.MVNOS.urban[[1]] <- predict(MVNOS.with.demographics[[1]], newdata = newdata.urban, type = "order")
pred.MVNOS.urban[[2]] <- predict(MVNOS.with.demographics[[2]], newdata = newdata.urban, type = "order")
pred.MVNOS.urban[[3]] <- predict(MVNOS.with.demographics[[3]], newdata = newdata.urban, type = "order")
dif.MVNOS.DID <- cbind(colMeans(rbind(t(apply(pred.MVNOS.urban[[1]]$o, c(2, 3), mean)), 
                                    t(apply(pred.MVNOS.urban[[2]]$o, c(2, 3), mean)), 
                                    t(apply(pred.MVNOS.urban[[3]]$o, c(2, 3), mean))) - 
                                rbind(t(apply(pred.MVNOS.rural[[1]]$o, c(2, 3), mean)), 
                                      t(apply(pred.MVNOS.rural[[2]]$o, c(2, 3), mean)), 
                                      t(apply(pred.MVNOS.rural[[3]]$o, c(2, 3), mean)))), 
                     HPDinterval(as.mcmc(rbind(t(apply(pred.MVNOS.urban[[1]]$o, c(2, 3), mean)), 
                                               t(apply(pred.MVNOS.urban[[2]]$o, c(2, 3), mean)), 
                                               t(apply(pred.MVNOS.urban[[3]]$o, c(2, 3), mean))) - 
                                           rbind(t(apply(pred.MVNOS.rural[[1]]$o, c(2, 3), mean)), 
                                                 t(apply(pred.MVNOS.rural[[2]]$o, c(2, 3), mean)), 
                                                 t(apply(pred.MVNOS.rural[[3]]$o, c(2, 3), mean))))))

# long-term partisanship
newdata.LDP <- newdata.DPJ <- newdata.JRP <- MVNOS.data
newdata.LDP[, 12:14] <- newdata.DPJ[, 12:14] <- newdata.JRP[, 12:14] <- 0
newdata.LDP$LDP <- 1
newdata.DPJ$DPJ <- 1
newdata.JRP$JRP <- 1
pred.MVNOS.LDP <- pred.MVNOS.DPJ <- pred.MVNOS.JRP <- list()
set.seed(5000)
pred.MVNOS.LDP[[1]] <- predict(MVNOS.with.demographics[[1]], newdata = newdata.LDP, type = "order")
pred.MVNOS.LDP[[2]] <- predict(MVNOS.with.demographics[[2]], newdata = newdata.LDP, type = "order")
pred.MVNOS.LDP[[3]] <- predict(MVNOS.with.demographics[[3]], newdata = newdata.LDP, type = "order")
pred.MVNOS.DPJ[[1]] <- predict(MVNOS.with.demographics[[1]], newdata = newdata.DPJ, type = "order")
pred.MVNOS.DPJ[[2]] <- predict(MVNOS.with.demographics[[2]], newdata = newdata.DPJ, type = "order")
pred.MVNOS.DPJ[[3]] <- predict(MVNOS.with.demographics[[3]], newdata = newdata.DPJ, type = "order")
pred.MVNOS.JRP[[1]] <- predict(MVNOS.with.demographics[[1]], newdata = newdata.JRP, type = "order")
pred.MVNOS.JRP[[2]] <- predict(MVNOS.with.demographics[[2]], newdata = newdata.JRP, type = "order")
pred.MVNOS.JRP[[3]] <- predict(MVNOS.with.demographics[[3]], newdata = newdata.JRP, type = "order")
dif.MVNOS.DPJ.LDP.2 <- cbind(colMeans(rbind(t(apply(pred.MVNOS.DPJ[[1]]$o, c(2, 3), mean)), 
                                          t(apply(pred.MVNOS.DPJ[[2]]$o, c(2, 3), mean)), 
                                          t(apply(pred.MVNOS.DPJ[[3]]$o, c(2, 3), mean))) - 
                                      rbind(t(apply(pred.MVNOS.LDP[[1]]$o, c(2, 3), mean)), 
                                            t(apply(pred.MVNOS.LDP[[2]]$o, c(2, 3), mean)), 
                                            t(apply(pred.MVNOS.LDP[[3]]$o, c(2, 3), mean)))), 
                           HPDinterval(as.mcmc(rbind(t(apply(pred.MVNOS.DPJ[[1]]$o, c(2, 3), mean)), 
                                                     t(apply(pred.MVNOS.DPJ[[2]]$o, c(2, 3), mean)), 
                                                     t(apply(pred.MVNOS.DPJ[[3]]$o, c(2, 3), mean))) - 
                                                 rbind(t(apply(pred.MVNOS.LDP[[1]]$o, c(2, 3), mean)), 
                                                       t(apply(pred.MVNOS.LDP[[2]]$o, c(2, 3), mean)), 
                                                       t(apply(pred.MVNOS.LDP[[3]]$o, c(2, 3), mean))))))
dif.MVNOS.JRP.LDP.2 <- cbind(colMeans(rbind(t(apply(pred.MVNOS.JRP[[1]]$o, c(2, 3), mean)), 
                                          t(apply(pred.MVNOS.JRP[[2]]$o, c(2, 3), mean)), 
                                          t(apply(pred.MVNOS.JRP[[3]]$o, c(2, 3), mean))) - 
                                      rbind(t(apply(pred.MVNOS.LDP[[1]]$o, c(2, 3), mean)), 
                                            t(apply(pred.MVNOS.LDP[[2]]$o, c(2, 3), mean)), 
                                            t(apply(pred.MVNOS.LDP[[3]]$o, c(2, 3), mean)))), 
                           HPDinterval(as.mcmc(rbind(t(apply(pred.MVNOS.JRP[[1]]$o, c(2, 3), mean)), 
                                                     t(apply(pred.MVNOS.JRP[[2]]$o, c(2, 3), mean)), 
                                                     t(apply(pred.MVNOS.JRP[[3]]$o, c(2, 3), mean))) - 
                                                 rbind(t(apply(pred.MVNOS.LDP[[1]]$o, c(2, 3), mean)), 
                                                       t(apply(pred.MVNOS.LDP[[2]]$o, c(2, 3), mean)), 
                                                       t(apply(pred.MVNOS.LDP[[3]]$o, c(2, 3), mean))))))
dif.MVNOS.JRP.DPJ.2 <- cbind(colMeans(rbind(t(apply(pred.MVNOS.JRP[[1]]$o, c(2, 3), mean)), 
                                          t(apply(pred.MVNOS.JRP[[2]]$o, c(2, 3), mean)), 
                                          t(apply(pred.MVNOS.JRP[[3]]$o, c(2, 3), mean))) - 
                                      rbind(t(apply(pred.MVNOS.DPJ[[1]]$o, c(2, 3), mean)), 
                                            t(apply(pred.MVNOS.DPJ[[2]]$o, c(2, 3), mean)), 
                                            t(apply(pred.MVNOS.DPJ[[3]]$o, c(2, 3), mean)))), 
                           HPDinterval(as.mcmc(rbind(t(apply(pred.MVNOS.JRP[[1]]$o, c(2, 3), mean)), 
                                                     t(apply(pred.MVNOS.JRP[[2]]$o, c(2, 3), mean)), 
                                                     t(apply(pred.MVNOS.JRP[[3]]$o, c(2, 3), mean))) - 
                                                 rbind(t(apply(pred.MVNOS.DPJ[[1]]$o, c(2, 3), mean)), 
                                                       t(apply(pred.MVNOS.DPJ[[2]]$o, c(2, 3), mean)), 
                                                       t(apply(pred.MVNOS.DPJ[[3]]$o, c(2, 3), mean))))))

## Figure A2 (simulation results on value ranking�fs relationship with long-term partisanship and demographic variables)
png("Figure_A2.png", width = 5, height = 6.8, units = "in", pointsize = 9, res = 1500)
layout(matrix(c(1:8), 4, 2, byrow = TRUE))
par(mar=c(4, 1, 3, 1))
dotchart(dif.MVNOS.female[c(1, 7:2), 1], labels = rev(value.names), pch = 20, xlim = c(-1, 1), 
         main = "Men \u2192 Women", xlab = "Change in Average Rank")
points(dif.MVNOS.female[c(1, 7:2), 1], c(1:7), pch = 19)
segments(dif.MVNOS.female[c(1, 7:2), 2], c(1:7), dif.MVNOS.female[c(1, 7:2), 3], c(1:7))
abline(v = 0, lty = 3)
dotchart(dif.MVNOS.age[c(1, 7:2), 1], labels = rev(value.names), pch = 20, xlim = c(-1, 1), 
         main = "Age 20s \u2192 70 and over", xlab = "Change in Average Rank")
points(dif.MVNOS.age[c(1, 7:2), 1], c(1:7), pch = 19)
segments(dif.MVNOS.age[c(1, 7:2), 2], c(1:7), dif.MVNOS.age[c(1, 7:2), 3], c(1:7))
abline(v = 0, lty = 3)
dotchart(dif.MVNOS.education[c(1, 7:2), 1], labels = rev(value.names), pch = 20, xlim = c(-1, 1), 
         main = "Junior high school \u2192 College", xlab = "Change in Average Rank")
points(dif.MVNOS.education[c(1, 7:2), 1], c(1:7), pch = 19)
segments(dif.MVNOS.education[c(1, 7:2), 2], c(1:7), dif.MVNOS.education[c(1, 7:2), 3], c(1:7))
abline(v = 0, lty = 3)
dotchart(dif.MVNOS.DID[c(1, 7:2), 1], labels = rev(value.names), pch = 20, xlim = c(-1, 1), 
         main = "DID population ratio 0 \u2192 1", xlab = "Change in Average Rank")
points(dif.MVNOS.DID[c(1, 7:2), 1], c(1:7), pch = 19)
segments(dif.MVNOS.DID[c(1, 7:2), 2], c(1:7), dif.MVNOS.DID[c(1, 7:2), 3], c(1:7))
abline(v = 0, lty = 3)
dotchart(dif.MVNOS.DPJ.LDP.2[c(1, 7:2), 1], labels = rev(value.names), pch = 20, xlim = c(-1, 1), 
         main = "LDP \u2192 DPJ", xlab = "Change in Average Rank")
points(dif.MVNOS.DPJ.LDP.2[c(1, 7:2), 1], c(1:7), pch = 19)
segments(dif.MVNOS.DPJ.LDP.2[c(1, 7:2), 2], c(1:7), dif.MVNOS.DPJ.LDP.2[c(1, 7:2), 3], c(1:7))
abline(v = 0, lty = 3)
dotchart(dif.MVNOS.JRP.LDP.2[c(1, 7:2), 1], labels = rev(value.names), pch = 20, xlim = c(-1, 1), 
         main = "LDP \u2192 JIP", xlab = "Change in Average Rank")
points(dif.MVNOS.JRP.LDP.2[c(1, 7:2), 1], c(1:7), pch = 19)
segments(dif.MVNOS.JRP.LDP.2[c(1, 7:2), 2], c(1:7), dif.MVNOS.JRP.LDP.2[c(1, 7:2), 3], c(1:7))
abline(v = 0, lty = 3)
dotchart(dif.MVNOS.JRP.DPJ.2[c(1, 7:2), 1], labels = rev(value.names), pch = 20, xlim = c(-1, 1), 
         main = "DPJ \u2192 JIP", xlab = "Change in Average Rank")
points(dif.MVNOS.JRP.DPJ.2[c(1, 7:2), 1], c(1:7), pch = 19)
segments(dif.MVNOS.JRP.DPJ.2[c(1, 7:2), 2], c(1:7), dif.MVNOS.JRP.DPJ.2[c(1, 7:2), 3], c(1:7))
abline(v = 0, lty = 3)
dev.off()

#### party-based comparison among Japanese political candidates ####
## Figure 6 (value structures and the distribution of value preferences among each parties�f 2016 HoC election candidates)
party.id.rearranged <- c(1, 2, 5, 3, 8, 4, 7, 6, 9, 10)
party.labels.UTAES <- c("LDP", "DP", "IfO", "CGP", "PJK", "JCP", "PLP&TYF", "SDP", "HRP", "Joint Candidates of\nOpposition Parties")
png("Figure_6.png", width = 6.4, height = 4.8, units = "in", pointsize = 9, res = 1500)
layout(matrix(c(1:12), 3, 4, byrow = TRUE))
par(mar=c(0, 0, 2, 0))
for (i in 1:10) {
  plot(density.circular(coord2rad(UTAES.obs.coords[UTAES.complete$party == party.id.rearranged[i], 1], 
                                  UTAES.obs.coords[UTAES.complete$party == party.id.rearranged[i], 2]), 
                        kernel = "vonmises", bw = 25), 
       main = party.labels.UTAES[i], axes = FALSE, shrink = 1.3, 
       xlim = c(-1.3, 1.3), ylim = c(-1.15, 1.45))
  points(UTAES.value.coords[, 1], UTAES.value.coords[, 2])
  points(0, 0, pch = 20, cex = 0.5)
  arrows.circular(coord2rad(mean(UTAES.obs.coords[UTAES.complete$party == party.id.rearranged[i], 1]), 
                            mean(UTAES.obs.coords[UTAES.complete$party == party.id.rearranged[i], 2])), 
                  y = (mean(UTAES.obs.coords[UTAES.complete$party == party.id.rearranged[i], 1]) ^ 2 + 
                         mean(UTAES.obs.coords[UTAES.complete$party == party.id.rearranged[i], 2]) ^ 2) ^ 0.5, 
                  length = 0.03)
  text(1.3, -1.15, labels = paste0("N = ", sum(UTAES.complete$party == party.id.rearranged[i])))
}
dev.off()

## Figure A3 (value structures and the distribution of value preferences of each party�fs winning candidates in the 2016 Japanese HoC election)
party.id.rearranged.winners <- c(1, 2, 5, 3, 4, 10)
party.labels.UTAES.winners <- c("LDP", "DP", "IfO", "CGP", "JCP", "Joint Candidates of\nOpposition Parties")
png("Figure_A3.png", width = 4.8, height = 3.2, units = "in", pointsize = 9, res = 1500)
layout(matrix(c(1:6), 2, 3, byrow = TRUE))
par(mar=c(0.5, 1, 2, 1))
for (i in 1:6) {
  plot(density.circular(coord2rad(UTAES.obs.coords[UTAES.complete$party == party.id.rearranged.winners[i] & 
                                                     UTAES.complete$RESULT == 1, 1], 
                                  UTAES.obs.coords[UTAES.complete$party == party.id.rearranged.winners[i] & 
                                                     UTAES.complete$RESULT == 1, 2]), 
                        kernel = "vonmises", bw = 25), 
       main = party.labels.UTAES.winners[i], axes = FALSE, shrink = 1.4, 
       xlim = c(-1.3, 1.3), ylim = c(-1.15, 1.45))
  points(UTAES.value.coords[, 1], UTAES.value.coords[, 2])
  points(0, 0, pch = 20, cex = 0.5)
  arrows.circular(coord2rad(mean(UTAES.obs.coords[UTAES.complete$party == party.id.rearranged.winners[i] & 
                                                    UTAES.complete$RESULT == 1, 1]), 
                            mean(UTAES.obs.coords[UTAES.complete$party == party.id.rearranged.winners[i] & 
                                                    UTAES.complete$RESULT == 1, 2])), 
                  y = (mean(UTAES.obs.coords[UTAES.complete$party == party.id.rearranged.winners[i] & 
                                               UTAES.complete$RESULT == 1, 1]) ^ 2 + 
                         mean(UTAES.obs.coords[UTAES.complete$party == party.id.rearranged.winners[i] & 
                                                 UTAES.complete$RESULT == 1, 2]) ^ 2) ^ 0.5, 
                  length = 0.03)
  text(1.3, -1.15, labels = paste0("N = ", sum(UTAES.complete$party == party.id.rearranged.winners[i] & 
                                                UTAES.complete$RESULT == 1)))
}
dev.off()

#### comparison between Japanese and U.S. voters ####
## loading Jacoby's (2014b) replication data
values.US <- read.table("ranked values, with names.txt", header = TRUE)  # value raking data
US.values <- values.US[, -1]
US.values.std <- t(apply(US.values, 1, scale))  # standadized value ranking for the MDPREF
info.US <- read.dta("info.dta")  # data containing respondents' party identification
merged.US <- merge(values.US, info.US, by = "caseid", all.x = TRUE)
merged.US$PID <- 3  # independent
merged.US$PID[merged.US$postpid7 >= 1 & merged.US$postpid7 <= 3] <- 1  # Democrat
merged.US$PID[merged.US$postpid7 >= 5 & merged.US$postpid7 <= 7] <- 2  # Republican

nrow(values.US)  # number of observations

## average ranking of values in the United States
ranking.mean.US <- colMeans(7 - US.values)
no1.share.US <- apply(US.values, 2, function(x) mean(x == 6)) * 100

## MDPREF for the U.S. data
mdpref.US <- ALSOS.MDPREF(data = US.values.std, tolerance = 0.01, max.iter = 25)
round(mdpref.US$fitvector, 3)
raw.freedom <- as.numeric(coord2rad(mdpref.US$decomp$v[1, 1], mdpref.US$decomp$v[1, 2]))
US.rotation.matrix <- matrix(c(cos(raw.freedom - pi / 2), sin(raw.freedom - pi / 2), 
                               -sin(raw.freedom - pi / 2), cos(raw.freedom - pi / 2)), 2, 2)
US.value.coords <- mdpref.US$decomp$v[, 1:2] %*% US.rotation.matrix
rownames(US.value.coords) <- value.names
US.obs.coords <- mdpref.US$subj.coord %*% US.rotation.matrix

## loading online survey data
online <- read.csv("online_survey_data.csv")
nrow(online)  # number of observations

## adjust the sample distribution of demographic variables
online.gender <- online$SEX - 1
online.age <- ifelse(online$AGE < 25, 1, 
                     ifelse(online$AGE < 30, 2, 
                            ifelse(online$AGE < 35, 3, 
                                   ifelse(online$AGE < 40, 4, 
                                          ifelse(online$AGE < 45, 5, 
                                                 ifelse(online$AGE < 50, 6, 
                                                        ifelse(online$AGE < 55, 7, 
                                                               ifelse(online$AGE < 60, 8, 
                                                                      ifelse(online$AGE < 65, 9, 
                                                                             ifelse(online$AGE < 70, 10, 
                                                                                    ifelse(online$AGE < 75, 11, 12)))))))))))
online.education <- ifelse(online$Q15 < 3, online$Q15, 
                           ifelse(online$Q15 < 6, 3, 4))
online.region <- ifelse(online$CHIIKI < 8, 1, 
                        ifelse(online$CHIIKI < 15, 2,  
                               ifelse(online$CHIIKI < 25, 3, 
                                      ifelse(online$CHIIKI < 31, 4, 
                                             ifelse(online$CHIIKI < 36, 5, 
                                                    ifelse(online$CHIIKI < 40, 6, 7))))))

# population distribution from the 2010 census
population.gender <- c(61841738, 65253007) / sum(c(61841738, 65253007))
population.age <- c(8369232, 6409612, 7290878, 8316157, 9732218, 8662804, 
                      7930296, 7515246, 8455010, 9643867, 7695811, 6276856) / 
  sum(c(8369232, 6409612, 7290878, 8316157, 9732218, 8662804, 
        7930296, 7515246, 8455010, 9643867, 7695811, 6276856))
population.education <- c(16756162, 41400268, 13187048, 17716535) / 
  sum(c(16756162, 41400268, 13187048, 17716535))
population.region <- c(14267, 43132, 23223, 20681, 7406, 2347, 14405) / 
  sum(c(14267, 43132, 23223, 20681, 7406, 2347, 14405))

# entropy balancing
ebalance.result <- ebalance(c(rep(0, nrow(online)), rep(1, nrow(online))), 
                            rbind(cbind(online.gender, diag(12)[online.age, ][, -1],
                                        diag(4)[online.education, ][, -1], 
                                        diag(7)[online.region, ][, -1]),
                                  tcrossprod(rep(1, nrow(online)),
                                             c(population.gender[-1], population.age[-1], 
                                               population.education[-1], population.region[-1]))))
online.weight <- ebalance.result$w

## average ranking of values in the online survey sample
online.raw.values <- online[, 2:8]
ranking.mean.online <- apply(online.raw.values, 2, weighted.mean, w = online.weight)
no1.share.online <- apply(online.raw.values, 2, function(x) weighted.mean(x == 1, w = online.weight)) * 100

## Figure A4 (average value ranks and the proportion choosing each value as their most important)
png("Figure_A4.png", width = 5.5, height = 4.2, units = "in", pointsize = 8, res = 1500)
layout(matrix(c(1:4), 2, 2, byrow = TRUE))
par(mar=c(4, 4, 3, 1))
plot(ranking.mean.UTAVS, no1.share.UTAVS, xlim = c(2, 6), ylim = c(0, 55), 
     main = "UTAVS (2014)", 
     xlab = "Average Ranks", ylab = "Chosen as Most Important (%)")
text(ranking.mean.UTAVS[1], no1.share.UTAVS[1], "Freedom", pos = 3)
text(ranking.mean.UTAVS[2], no1.share.UTAVS[2], "Equality", pos = 3)
text(ranking.mean.UTAVS[3], no1.share.UTAVS[3], "Economic\nstability", pos = 3)
text(ranking.mean.UTAVS[4], no1.share.UTAVS[4], "Morality", pos = 3)
text(ranking.mean.UTAVS[5], no1.share.UTAVS[5], "Self-reliance", pos = 3)
text(ranking.mean.UTAVS[6], no1.share.UTAVS[6], "Social\norder", pos = 3)
text(ranking.mean.UTAVS[7] - 0.1, no1.share.UTAVS[7], "Patriotism", pos = 1)
plot(ranking.mean.online, no1.share.online, xlim = c(2, 6), ylim = c(0, 55), 
     main = "Online Survey (2017)", 
     xlab = "Average Ranks", ylab = "Chosen as Most Important (%)")
text(ranking.mean.online[1], no1.share.online[1], "Freedom", pos = 3)
text(ranking.mean.online[2], no1.share.online[2], "Equality", pos = 1)
text(ranking.mean.online[3] - 0.2, no1.share.online[3], "Economic\nstability", pos = 1)
text(ranking.mean.online[4], no1.share.online[4], "Morality", pos = 4)
text(ranking.mean.online[5], no1.share.online[5], "Self-reliance", pos = 3)
text(ranking.mean.online[6], no1.share.online[6], "Social\norder", pos = 3)
text(ranking.mean.online[7], no1.share.online[7], "Patriotism", pos = 2)
plot(ranking.mean.US, no1.share.US, xlim = c(2, 6), ylim = c(0, 50), 
     main = "2006 CCES", 
     xlab = "Average Ranks", ylab = "Chosen as Most Important (%)")
text(ranking.mean.US[1], no1.share.US[1], "Freedom", pos = 3)
text(ranking.mean.US[2], no1.share.US[2], "Equality", pos = 1)
text(ranking.mean.US[3] - 0.2, no1.share.US[3], "Economic\nsecurity", pos = 3)
text(ranking.mean.US[4], no1.share.US[4], "Morality", pos = 4)
text(ranking.mean.US[5], no1.share.US[5], "Individualism", pos = 3)
text(ranking.mean.US[6] - 0.2, no1.share.US[6], "Social\norder", pos = 1)
text(ranking.mean.US[7], no1.share.US[7], "Patriotism", pos = 1)
dev.off()

## Figure A5 (value structures and the distribution of value preferences)
png("Figure_A5.png", width = 6.4, height = 3.2, units = "in", pointsize = 9, res = 1500)
layout(matrix(c(1:2), 1, 2, byrow = TRUE))
par(mar=c(0.5, 0.5, 3, 0.5))
set.seed(12345)
plot(density.circular(coord2rad(UTAVS.obs.coords[, 1], UTAVS.obs.coords[, 2]), 
                      kernel = "vonmises", bw = 25), main = "Japanese Voters (2014)", axes = FALSE, shrink = 1.2)
points(UTAVS.value.coords[, 1], UTAVS.value.coords[, 2])
text(UTAVS.value.coords[, 1], UTAVS.value.coords[, 2], labels = value.names, cex = 0.8, pos = 3)
arrows.circular(coord2rad(mean(UTAVS.obs.coords[, 1]), mean(UTAVS.obs.coords[, 2])), 
                y = (mean(UTAVS.obs.coords[, 1]) ^ 2 + mean(UTAVS.obs.coords[, 2]) ^ 2) ^ 0.5, 
                length = 0.05)
sampled.id <- sample(c(1:nrow(UTAVS.complete)), 300)
ticks.circular(coord2rad(UTAVS.obs.coords[, 1], UTAVS.obs.coords[, 2])[sampled.id])
text(1.2, -1.2, labels = paste0("N = ", nrow(UTAVS.obs.coords)))
plot(density.circular(coord2rad(-US.obs.coords[, 1], US.obs.coords[, 2]), 
                      kernel = "vonmises", bw = 25), main = "U.S. Voters (2006)", axes = FALSE, shrink = 1.2)
points(-US.value.coords[, 1], US.value.coords[, 2])
text(-US.value.coords[, 1], US.value.coords[, 2], 
     labels = c("Freedom", "Equality", "Economic security", "Morality", 
                "Individualism", "Social order", "Patriotism"), 
     cex = 0.8, pos = c(3, 3, 3, 3, 2, 3, 3))
arrows.circular(coord2rad(mean(-US.obs.coords[, 1]), mean(US.obs.coords[, 2])), 
                y = (mean(-US.obs.coords[, 1]) ^ 2 + mean(US.obs.coords[, 2]) ^ 2) ^ 0.5, 
                length = 0.05)
sampled.id <- sample(c(1:nrow(values.US)), 300)
ticks.circular(coord2rad(-US.obs.coords[, 1], US.obs.coords[, 2])[sampled.id])
text(1.2, -1.2, labels = paste0("N = ", nrow(US.obs.coords)))
dev.off()

## Figure A6 (value structures and the distribution of value preferences of U.S. voters by party identification)
party.labels.US <- c("Democrats", "Republicans", "Independents")
png("Figure_A6.png", width = 4.8, height = 1.6, units = "in", pointsize = 9, res = 1500)
layout(matrix(c(1:3), 1, 3))
par(mar=c(0.5, 1, 2, 1))
for (i in 1:3) {
  plot(density.circular(coord2rad(-US.obs.coords[merged.US$PID == i, 1], 
                                  US.obs.coords[merged.US$PID == i, 2]), 
                        kernel = "vonmises", bw = 25), 
       main = party.labels.US[i], axes = FALSE, shrink = 1.75)
  points(-US.value.coords[, 1], US.value.coords[, 2])
  points(0, 0, pch = 20, cex = 0.5)
  arrows.circular(coord2rad(mean(-US.obs.coords[merged.US$PID == i, 1]), 
                            mean(US.obs.coords[merged.US$PID == i, 2])), 
                  y = (mean(-US.obs.coords[merged.US$PID == i, 1]) ^ 2 + 
                         mean(US.obs.coords[merged.US$PID == i, 2]) ^ 2) ^ 0.5, 
                  length = 0.03)
  text(1.4, -1.4, labels = paste0("N = ", sum(merged.US$PID == i)))
}
dev.off()